PMF不是画出来就算数：从收敛、重加权到2D自由能面的物理判据

PMF不是画出来就算数：从收敛、重加权到2D自由能面的物理判据很多人第一次做 PMF 时，最容易掉进一个坑：图是画出来了，但物理上并不一定成立。问题在于，能画出来，和能不能当成平衡自由能解释，是两回事。这篇文章只回答几个更基础、也更容易出错的问题：已有数据什么时候足够支持 PMF，什么时候只能报局部结果，什么时候必须重加权，什么时候 2D 图虽然能画，但其实不该把它写成“收敛的自由能面”。结论 PMF 的定义本身并不难，真正困难的是采样是否真的支持这个定义。无偏 MD 确实可以直接给自由能，但前提是分析段已经平稳，而且目标坐标空间被充分访问；只要存在偏置、约束、umbrella 或多窗口合并，就不能跳过重加权。 2D PMF 不是“多画一个维度”那么简单，而是对采样混合提出了更高要求。如果某些区域从来没被访问过，任何后处理都不能把真实自由能补出来；因此，很多时候你真正能安全报告的，并不是全局 PMF，而是局部 PMF、条件分布或状态占据。 PMF 到底是什么对一个集合变量 $\xi$，平衡自由能剖面定义为： \[F(\xi) = -k_B T \ln P(\xi) + C\] 如果有两个集合变量 $\xi,\eta$，对应的二维自由能面就是： \[F(\xi,\eta) = -k_B T \ln P(\xi,\eta) + C\] 公式的通俗解释这两个式子真正表达的是一句很朴素的话：某个状态如果在平衡系综里更常出现，它的自由能就更低。所以，问题的核心从来不是“会不会取负对数”，而是你算出来的 $P(\xi)$ 或 $P(\xi,\eta)$ 到底是不是平衡分布，这个分布覆盖的是全局空间还是只覆盖了一个局部盆地，以及每个 bin 里到底有多少有效独立样本。这三件事，才真正决定了你的 PMF 能不能被当成物理结果来解释。在后面的例子里，我会经常用 P2 和 Z 这两个符号。这里可以先把它们通俗地理解成两类常见坐标：P2 代表某种取向序参量，也就是“分子更偏向平躺、倾斜还是竖直”的量化描述；Z 代表某种位置坐标，例如分子相对于界面、膜中心或参考平面的距离。你完全可以把它们替换成自己体系里真正关心的两个集合变量。什么叫“物理上正确”的 PMF 如果想让一条 PMF 在文章里站得住脚，至少要同时满足四件事：数据来自同一个目标系综用来分析的轨迹段已经进入平稳区你关心的坐标范围内发生了足够的往返跃迁误差估计使用的是有效样本数，不是总帧数只要这四条里缺一条，图可能仍然能画出来，但解释时就必须明显降级。第一关：是不是同一个统计系综这一点最容易被忽视。如果所有数据都来自同一统计系综，也就是温度一致、压力设置一致、力场和拓扑一致、体系组成与边界条件一致，同时没有额外偏置或约束，那么这些轨迹才有资格被当作同一个平衡分布的样本来合并分析。那么你可以直接从直方图或核密度估计（KDE）得到 $P(\xi)$，再转成自由能。但只要出现下面任一种情况，就不能把所有帧直接混在一起做直方图：情况为什么不能直接混合对某个坐标加了 umbrella 势采样分布已经被显式改权，不再对应原始无偏分布加了位置约束或取向约束体系访问相空间的方式被限制，直方图不再代表自然占据做过 steered MD 或 pulling 轨迹带有外场驱动，不能直接当成平衡样本合并了不同温度的数据不同温度对应不同平衡分布，不能简单拼接合并了不同哈密顿量或不同参数的数据势能面本身不同，统计权重自然也不同这时你要处理的已经不是“无偏概率”，而是“被改权重后的采样概率”。必须重加权，常见工具就是 WHAM、MBAR，或者更一般的重加权流程。第二关：轨迹是不是已经进入平稳区很多 PMF 最大的问题，不是采样短，而是前半段根本还没平衡。比如系统一开始从某个强行构建的初始构型出发，前几十纳秒甚至更久都还在弛豫。如果把这一段直接并进统计，得到的就不是平衡分布，而是“初始条件残留 + 平衡波动”的混合物。一个实用做法，是先做平衡段检测，再决定从哪里开始统计。常用工具是 pymbar.timeseries。这里输入的数据，不是什么特殊格式文件，而是某个集合变量随时间变化的一列数据，最常见的就是 P2(t) 或 Z(t) 这样的时间序列： python - <<'PY' from pymbar import timeseries import numpy as np P2_t = np.loadtxt('P2_t.dat') t0, g, Neff = timeseries.detect_equilibration(P2_t, nskip=10) print(t0, g, Neff) PY 如果你手里保存的是多列文件，例如同一份文件里同时有时间、P2 和 Z，那就应该先把你想分析的那一列取出来，再送进 detect_equilibration()，而不是把整张表不加区分地直接读进去。这里最值得报告的，不是“我跑了多少 ns”，而是平衡起点 $t_0$、统计低效因子 $g$ 和有效样本数 $N_{\mathrm{eff}}$。真正决定误差条大小的，是独立样本有多少，不是帧有多少。很多时候看起来“已经有几十万帧”，但如果自相关很强，真正能用于统计判断的独立样本可能并不多。第三关：有没有真正发生“来回走动” 这是判断 PMF 是否可信的核心。真正有用的判断，不是“分布看起来挺宽”，而是体系有没有在你关心的几个主要状态之间真正来回走动，也就是是否发生了足够多的往返跃迁（round trips）。对 1D 和 2D PMF，要求到底差在哪里目标至少要看到什么不能轻易下的结论 1D PMF 主要盆地被多次访问，盆地之间有往返跃迁，不同重复给出相近边缘分布只有单盆地波动时，不应宣称得到全局 PMF 2D PMF 两个坐标都被实质性访问，且在固定第一维时第二维也能混合，不同区域之间整体连通如果第二维几乎没动，或固定某一维后另一维几乎不跨峰，就不应宣称得到全局 2D 自由能面如果体系只在一个盆地附近晃动，那么你当然也能画出一条曲线，但那更接近“局部热涨落的自由能近似”，而不是全局 PMF。二维情况则更严格，因为它要求你不仅采到 $\xi$，还要在不同 $\eta$ 条件下把 $\xi$ 也采匀；一旦第二维只是窄范围波动，这张 2D 图通常就只能算局部地形。一个最常见的误区：能画 2D，不等于应该发 2D 很多人会这样做：选两个坐标，做二维直方图，再对联合概率取负对数，最后得到一张彩色图。从程序角度看完全没问题，但从物理角度看，可能只说明一件事：你的轨迹在一个局部区域里留下了很多点。这时真正应该问的，不是“图是不是好看”，而是三个更扎实的问题。第一，第二维是不是只覆盖了一个很窄的范围；如果是，那么 2D 图只是把局部波动展开成二维，并没有真正回答更大的自由能问题。第二，高自由能区域是“真的高”，还是“根本没采到”；没有访问到的格点，在视觉上很容易被误读成高能区，但统计学上它可能只是空白区。第三，盆地之间的通道是物理能垒，还是统计断裂；如果两个盆地中间几乎没有过渡点，你看到的未必是高能屏障，也可能只是采样没有连通，更专业地说，就是这些区域之间缺少足够的统计连通性。如果这些问题答不上来，最稳妥的表述通常不是“得到了全局 2D PMF”，而是把口径主动降到“局部 2D 自由能地形”“条件分布 $P(\xi\mid\eta)$”或者“已结合区间内的取向自由能”。什么时候无偏 MD 足够无偏 MD 适合回答的问题，其实比很多人想象得更有限，但也更扎实。与其笼统地说“能不能算 PMF”，不如先区分你到底想回答哪一类问题。目标无偏 MD 的适用性更合适的表述单个坐标的 1D 边缘自由能较好 1D PMF 某个局部区域内的自由能起伏较好局部 PMF 分箱后的状态占据比较较好条件分布或占据统计跨多个盆地的全局自由能谨慎只有在多次跨盆地跃迁后才可报告同时含位置与取向的 2D 自由能面很谨慎通常先降级为局部 2D 或条件分布含解离、再结合、重排等慢过程很谨慎往往需要增强采样支撑如果你的无偏轨迹从头到尾都没有离开某个状态盆地，那么最合理的结论不是“体系没有别的态”，而是：当前采样没有能力回答这个问题。什么时候必须用 WHAM 或 MBAR 这个判断其实很干脆：只要采样权重被改过，就要重加权。与其把这一条说成一句口号，不如直接看常见场景：场景能不能直接做直方图推荐处理同一无偏 MD 可以直方图或 KDE umbrella 窗口不可以 WHAM 或 MBAR 多温度数据合并不可以 MBAR 有约束或 pulling 不可以显式重加权多个偏置窗口做 2D 分布不可以先去偏，再做联合分布如果你手上已有沿某个坐标布置好的 umbrella 窗口，那么它们通常足够支持可靠的 1D PMF。至于能不能进一步得到 2D PMF，要看另一个坐标在每个窗口里是不是也混合得足够好。主坐标被偏置采到，并不自动意味着旁观变量也已经收敛，这一点在实际分析里经常被误判。一个非常实用的判断：你到底能安全声称什么诊断结果最稳妥的说法只有一个局部盆地被采到局部自由能或局部涨落 1D 有多次跨峰跃迁，重复一致可以报告 1D PMF 2D 中第二维很窄只报告条件分布或局部 2D 地形 umbrella 在主坐标重叠良好，但副坐标混合差主坐标 PMF 可信，2D 结果仅作定性参考每个窗口内副坐标多次跨峰，重复一致可以认真讨论 2D PMF 这张表背后的原则其实很简单：结论的口径，必须和采样能力匹配。很多结果并不是“完全不能发”，而是应该主动把口径降到“局部 PMF”“条件分布”或者“占据统计”这一层，这样反而更稳。收敛不能只看“曲线变平” 很多人判断收敛时，只看 PMF 曲线后半段是不是“不怎么变了”。这远远不够，因为一条表面平滑的曲线，可能只是建立在高度相关、重复不一致、或者根本没有跨盆地跃迁的数据上。更可靠的收敛证据链更可靠的判断，通常要把下面几类证据合在一起看：先看结果会不会随时间继续漂，也就是是否仍在发生系统性漂移；再看不同重复是否支持同一组物理结论；接着看你到底有多少真正独立的样本；最后再确认主要状态之间有没有真正发生来回切换，也就是是否存在足够的往返跃迁。时间分块分析：把前 1/3、前 2/3 和全部数据分别算一次 PMF。这样做的目的，不是为了多画几条线，而是看结果会不会继续变。如果主要盆地位置、相对深度和势垒高度还在系统性漂移，那就说明体系还在持续演化、尚未真正稳定下来，此时“看起来平滑”并不等于已经收敛。重复一致性：不同重复轨迹给出的分布或 PMF 应该大体一致。这里最重要的不是三条线能不能完全重合，而是它们是否支持同一个物理结论。如果不同重复之间差异明显，最常见的解释不是“体系本来就这样”，而是混合仍然不足，也就是每条轨迹还在各自记着不同的初始路径。自相关分析：报告 $g$ 和 $N_{\mathrm{eff}}$，确认自己不是在用几十万帧去假装拥有几十万个独立样本。连续轨迹里的相邻帧往往很像，所以“帧数很多”不等于“信息很多”。这一步本质上是在修正相关样本导致的误差低估，也就是给误差条去水分，说明到底有多少真正能独立贡献统计信息的数据点。跃迁计数：主要盆地之间要有实质性的往返，而不是只在一个盆地里高频抖动。很多人看到时间序列很活跃，就以为体系采样得很好，但如果这些波动始终发生在同一个局部盆地里，那么关键状态之间的相对自由能差其实还没有被真正比较过。没有跨盆地跃迁时，很多相对自由能差并不稳。窗口重叠：对 umbrella 来说，相邻窗口必须足够连通。如果相邻窗口之间几乎没有共同覆盖的区域，WHAM 或 MBAR 就很难把整条 PMF 稳稳地拼起来。这时数学上虽然还能算，物理上却可能只是把几段彼此脱节的局部结果硬接在一起；更规范地说，就是窗口之间缺少足够的概率分布重叠。 umbrella 数据至少要看什么对于 umbrella，gmx wham 的常规检查项很重要： gmx wham -it tpr-files.dat -if pullf-files.dat -o pmf.xvg -hist hist.xvg -ac 这里至少要看三件事，而且最好把它们理解成“这条 PMF 能不能被顺畅接起来”的三个层次检查：相邻窗口直方图有没有足够重叠。这是最基础的一关。如果相邻窗口几乎不相交，那么后处理再漂亮，也只是把统计上彼此脱节的区间强行缝在一起，整条曲线会缺少真正的连接。自相关时间是不是已经大到接近单窗口长度。这一步是在问：单个窗口里到底有没有采到足够多的独立信息。如果一个窗口里有效独立样本本来就很少，那么它对整条 PMF 的贡献会既不稳定又很难估误差；此时窗口数量再多，也不等于每个窗口都真的达到局部统计稳定。不同窗口拼起来后有没有明显断链。所谓断链，不一定表现成肉眼可见的大跳跃，也可能表现为某些区间误差异常、重复不一致，或者对分析参数极其敏感。如果一条 PMF 只要稍微改一下 bin、平滑或截断方式就明显变样，那通常不是“图画风不同”，而是底层采样还不够扎实。如果某些窗口几乎没有重叠，或者窗口内采样时间和自相关时间是一个量级，那这套 PMF 就很难让人放心。 2D PMF 什么时候才值得做更关键的问题是：什么时候做 2D PMF 比做 1D 或条件分布更有信息增益。通常至少要同时满足三点：两个坐标都对应你真正关心的慢过程，这两个坐标在数据里都被实质性采样到了，而且在固定第一维时第二维不是“卡死”的，也就是没有被困在某个狭窄取值范围里。少了其中任何一条，二维分析带来的往往不是新信息，而是新噪声。如果不满足，2D 往往只会带来两个后果：图更花哨，误差更大。因为二维一上来就会遭遇“维数灾难”：格点数一多，平均到每个 bin 的有效样本数会迅速下降，空 bin 和噪声会明显增加。所以，在下面这些情况下，不做 2D 反而更专业：如果第二维只是辅助解释变量，如果第二维的采样范围很窄，如果第二维的混合时间明显比单窗口长度更长，或者你的核心结论本质上靠 1D 就已经成立，那么继续硬做 2D 往往只会增加图的复杂度，而不会提高结论的可信度。还有一个细节：有些序参量自带“几何熵” 如果你用的是角度、取向序参量，或者由角度变换得到的量，那么要小心一个问题：原始分布里可能混进了变量测度本身带来的偏置。最直观的例子就是方向相关变量。即使体系完全各向同性，某些取向序参量的概率分布也未必是均匀的。这意味着直接计算 \[F(\xi) = -k_B T \ln P(\xi) + C\] 得到的可能既包含真实相互作用偏好，也包含“随机几何本来就更容易落在某些值附近”的贡献。这时最常见的处理方式有两种：报告方式含义适合的讨论场景原始 PMF 包含变量测度带来的几何熵讨论状态占据、总体分布相对参考分布的超额自由能更突出相互作用导致的偏好讨论取向偏好、界面诱导效应这不是所有体系都必须做，但如果你的核心结论高度依赖“取向偏好”，那这个问题最好提前想清楚。否则读者看到的“最低谷”，有一部分可能只是变量定义自带的几何效应，而不全是体系相互作用本身。一个面向实战的工作流 graph TB A["拿到已有轨迹"] --> B["先分清：无偏数据还是有偏数据"] B --> C["确定目标：1D、局部2D、还是全局2D"] C --> D["检测平衡段：t0、g、Neff"] D --> E["检查跃迁、重复一致性、窗口重叠"] E --> F{"采样是否支持目标结论"} F -->|支持| G["报告 PMF，并给出误差与收敛证据"] F -->|部分支持| H["降级为局部 PMF、条件分布或状态占据"] F -->|不支持| I["补采样或重新设计增强采样方案"] 这个流程最重要的一步，不是“画图”，而是中间那个判断：采样能力到底支不支持你想说的话。真正成熟的分析，不是把所有图都画出来，而是知道哪些图值得认真解释，哪些图只能当辅助材料。结果该怎么讲，才更站得住脚一张自由能图要站得住脚，关键不在于修饰，而在于先把哪里可信、哪里还不能多说讲清楚：先说明平衡段和有效样本是怎么处理的。如果一开始就交代你已经剔除了前期非平衡部分，并且按相关性修正了有效样本数，读者会更容易接受后面的自由能结果，因为他知道这些曲线不是把所有帧不加区分地堆出来的。再说明 1D 结果为什么可信。如果主要状态之间已经出现多次往返跃迁，而且不同重复支持同一个结论，那么这时去讨论 1D PMF 的相对高低才更有底气，因为它背后有明确的动力学采样证据。谈到 2D 结果时主动限定范围。如果二维图只有一部分区域采样得比较扎实，那就只讨论那一部分，把它明确写成局部自由能地形或条件分布。这样做不会削弱文章，反而会让读者觉得你的判断更稳。对空白区和混合不足区保持克制。没有访问到的区域就不要硬解释，混合明显不足的方向也不要勉强下定量结论。这样做不是示弱，而是在保护结论的可信度。这种写法的价值不在于“更谨慎”，而在于把真正确定的部分讲扎实，把暂时不能确定的部分老老实实留白。最后总结 PMF 真正难的地方，从来不是软件命令，而是你是否对“这张图能回答什么问题”有清醒判断。无偏 MD 确实可以直接给自由能，但前提是轨迹分析段已经平稳、混合、可重复。如果连主要状态之间的往返都没有发生，那么图上看到的更多只是局部波动，而不是可以放心解释的全局自由能。只要数据里存在偏置、约束、umbrella 或多窗口拼接，就必须认真做重加权。这不是后处理里的可选美化步骤，而是把“被改过权重的采样”还原成目标分布所必需的物理操作。 2D PMF 的门槛显著高于 1D PMF，因为它要求两个坐标都被充分访问，而且在固定其中一维时另一维也要发生足够混合。很多 1D 看起来已经稳定的数据，一到二维分析就会暴露出空白区、断裂区和高噪声问题。没采到就是没采到，后处理不能替代真实采样。无论是更平滑的直方图、更复杂的重加权，还是更漂亮的二维彩图，都不能凭空恢复从未被访问过的状态或通道。当采样只支持局部结论时，老老实实报告局部结论，反而更有说服力。把结果写成局部 PMF、条件分布或状态占据，通常比强行宣称“全局自由能面已经收敛”更专业，也更经得起追问。如果把这套判断标准先建立起来，你之后无论做无偏 MD、umbrella、metadynamics，还是更复杂的多维自由能分析，很多技术决策都会清楚得多。

Molecular Dynamics · 2026-03-31

靶向分子动力学（TMD）：用RMSD约束引导蛋白质构象转变

Molecular Dynamics · 2025-10-11

自适应偏置力（ABF）方法详解

自适应偏置力（ABF）方法详解一、ABF方法的基本原理自适应偏置力（Adaptive Biasing Force, ABF）是一种用于计算自由能曲面（PMF）的增强采样方法。它的核心思想是：通过实时计算并施加一个抵消系统平均力的偏置力，使分子能够在反应坐标上自由扩散，从而加速采样。基本方程对于一个集合变量（collective variable, CV）$\xi$，系统在 $\xi$ 方向上受到的瞬时力为 $F(\xi)$。ABF方法通过累积统计，估算出在 $\xi$ 处的平均力 $\langle F(\xi) \rangle$： \[\langle F(\xi) \rangle = -\frac{\mathrm{d}A(\xi)}{\mathrm{d}\xi}\] 其中 $A(\xi)$ 是沿着 $\xi$ 的自由能（PMF）。 ABF的策略：在模拟过程中，实时施加一个偏置力 $F_{bias}(\xi) = -\langle F(\xi) \rangle$，使得分子在 $\xi$ 方向上受到的净力接近零，从而能够自由地在整个 $\xi$ 范围内扩散。瞬时力的计算：从原子力到集合变量的投影关键问题：MD引擎（如NAMD、GROMACS）计算的是原子间的相互作用力 $\mathbf{F}_i$（作用在每个原子 $i$ 上），但ABF需要的是沿着集合变量 $\xi$ 的广义力 $F(\xi)$。如何将原子力转换为CV方向的力？答案：通过链式法则投影。集合变量 $\xi$ 通常是原子坐标 ${\mathbf{r}_i}$ 的函数，即 $\xi = \xi(\mathbf{r}_1, \mathbf{r}_2, \ldots, \mathbf{r}_N)$。瞬时力通过以下公式计算： \[F(\xi) = -\sum_{i=1}^{N} \mathbf{F}_i \cdot \frac{\partial \xi}{\partial \mathbf{r}_i}\] 物理意义： $\frac{\partial \xi}{\partial \mathbf{r}_i}$ 是CV对第 $i$ 个原子坐标的梯度，表示该原子沿哪个方向运动会增加 $\xi$ 的值 $\mathbf{F}_i \cdot \frac{\partial \xi}{\partial \mathbf{r}_i}$ 是原子 $i$ 受到的力在CV方向上的投影分量负号是因为力的定义（$\mathbf{F} = -\nabla U$）具体例子：在本文中，CV是小分子沿膜法线（z轴）的位置，即 $\xi = z_{molecule}$。此时： $\frac{\partial \xi}{\partial \mathbf{r}_i} = (0, 0, 1)$ 只有z分量非零 $F(\xi) = -F_{i,z}$ 只需提取分子受力的z分量实际实现：每个MD时间步，MD引擎计算所有原子受到的力 ${\mathbf{F}_i}$ Colvars模块（NAMD）或相应的插件（GROMACS）实时计算：当前的CV值 $\xi(t)$ CV的梯度 ${\partial\xi/\partial\mathbf{r}_i}$ 瞬时广义力 $F(\xi,t)$ 累积到直方图：将 $F(\xi,t)$ 加到对应 $\xi$ 网格点的累积和中计算平均力：$\langle F(\xi) \rangle = \frac{1}{N_{samples}(\xi)} \sum_{t:\xi(t)\approx\xi} F(\xi,t)$ 施加偏置：在下一个时间步，对相关原子施加偏置力 $\mathbf{F}_{bias,i} = -\langle F(\xi) \rangle \cdot \frac{\partial \xi}{\partial \mathbf{r}_i}$ 技术细节： ABF使用分层网格将CV空间离散化（如每0.01 nm一个网格点）为避免初期统计不准确，通常设置最小采样阈值（如每个网格点至少100次访问）才开始施加偏置力偏置力的施加使用渐进式缩放（ramp），从0逐渐增加到1，避免非平衡效应自由能的恢复模拟结束后，通过对累积的平均力进行积分，即可恢复自由能曲面： \[A(\xi) = A(\xi_0) - \int_{\xi_0}^{\xi} \langle F(\xi') \rangle \mathrm{d}\xi'\] 二、ABF的窗口策略与边界处理为什么需要分窗口？虽然理论上ABF可以在整个反应坐标范围内一次性进行（全局ABF），但在实际应用中，当自由能曲面存在高能垒时，全局ABF会遇到严重的采样问题：能垒区域采样不足：分子很难跨越高能垒区域，导致这些区域的平均力估计不准确收敛极慢：即使施加了偏置力，分子在能垒区域的停留时间仍然很短，需要极长的模拟时间才能充分采样解决方案：将整个反应坐标范围划分为多个重叠的窗口（stratification），在每个窗口内独立进行ABF采样，最后将各窗口的PMF拼接起来。窗口的定义每个窗口由以下参数定义：窗口范围 $[\xi_{min}, \xi_{max}]$：CV允许的取值范围窗口宽度：$\Delta\xi = \xi_{max} - \xi_{min}$（本文中为0.4 nm）窗口中心：$\xi_{center} = (\xi_{min} + \xi_{max})/2$ 相邻窗口的间隔：中心点之间的距离（本文中为0.1 nm）例如，在本文中：窗口1：$[-0.2, +0.2]$ nm，中心在 0 nm 窗口2：$[-0.1, +0.3]$ nm，中心在 +0.1 nm 窗口3：$[0.0, +0.4]$ nm，中心在 +0.2 nm … 边界的处理方式 ABF方法对窗口边界的处理与umbrella sampling有本质区别： 1. 无强制约束的边界 ABF不在窗口边界施加强制约束势。当CV的值 $\xi$ 处于窗口范围 $[\xi_{min}, \xi_{max}]$ 内时：正常施加偏置力：$F_{bias}(\xi) = -\langle F(\xi) \rangle$ 正常采样和累积统计：该位置的构象被记录用于平均力的估算当 $\xi$ 超出窗口范围时：停止施加偏置力：不再对系统施加ABF偏置停止采样：该位置的构象不被记录模拟继续运行：系统仍然正常演化，只是不参与当前窗口的统计 2. 可选的软约束势（wall potential）为了防止分子过度偏离窗口范围，可以在边界外侧添加一个软约束势（也称为wall potential或restraining potential）： \[U_{wall}(\xi) = \begin{cases} \frac{k}{2}(\xi - \xi_{max})^2 & \text{if } \xi > \xi_{max} + \delta \\ 0 & \text{if } \xi_{min} - \delta \leq \xi \leq \xi_{max} + \delta \\ \frac{k}{2}(\xi - \xi_{min})^2 & \text{if } \xi < \xi_{min} - \delta \end{cases}\] 其中： $k$ 是弹簧常数（通常为10-100 kcal/mol/Å²） $\delta$ 是缓冲区宽度（通常至少为一个网格间距）关键特点：约束势的作用范围应比窗口范围更宽（$\delta > 0$），确保在窗口边界处没有突变约束势是柔和的（软约束），不会强制将分子”锁死”在某个位置与Umbrella Sampling的对比特性 ABF Umbrella Sampling 窗口定义定义边界范围 $[\xi_{min}, \xi_{max}]$ 定义中心点 $\xi_0$ 约束方式无强制约束（或软约束）强制谐振子势 $\frac{k}{2}(\xi-\xi_0)^2$ 分子运动在整个窗口内自由扩散被”拴”在中心点附近，受弹簧限制偏置力动态调整，实时抵消平均力静态谐振子势后处理不需要，直接积分平均力得PMF 需要WHAM等方法去除偏置先验知识不需要知道自由能形状需要预估PMF形状来设置弹簧常数窗口重叠不强制要求（但推荐）必须重叠，否则WHAM无法拼接三、窗口的拼接与PMF的构建重叠区域的作用虽然ABF在理论上不强制要求窗口重叠（因为平均力是连续的），但在实践中高度推荐使用重叠窗口，原因如下：提高统计精度：重叠区域被两个窗口同时采样，提供了交叉验证平滑过渡：减少拼接时的不连续性检测采样质量：如果两个窗口在重叠区域的PMF差异很大，说明采样不充分拼接算法详解 ABF窗口拼接的核心挑战在于：每个窗口独立模拟得到的PMF只是相对值（积分常数未定），需要通过重叠区域将它们”对齐”到同一个能量基准上。步骤1：对每个窗口内的平均力进行积分对于第 $i$ 个窗口（范围 $[\xi_i^{min}, \xi_i^{max}]$），从下边界开始积分平均力： \[A_i(\xi) = -\int_{\xi_i^{min}}^{\xi} \langle F_i(\xi') \rangle \mathrm{d}\xi', \quad \xi \in [\xi_i^{min}, \xi_i^{max}]\] 注意：这里人为设定 $A_i(\xi_i^{min}) = 0$，所以 $A_i(\xi)$ 只是窗口内的相对PMF 积分通常使用数值方法（如梯形法则或辛普森法则）如果平均力在某些点采样不足，可能需要平滑处理（如样条插值）步骤2：在重叠区域对齐相邻窗口对于相邻的窗口 $i$ 和 $i+1$，它们的重叠区域是 $[\xi_{i+1}^{min}, \xi_i^{max}]$。在这个区域内，两个窗口都提供了PMF估计：$A_i(\xi)$ 和 $A_{i+1}(\xi)$。目标：找到一个偏移常数 $\Delta A_i$，使得 $A_i(\xi) + \Delta A_i \approx A_{i+1}(\xi)$ 在重叠区域内尽可能一致。方法1：简单平均法 $\Delta A_i = \frac{1}{N_{overlap}} \sum_{\xi \in overlap} [A_{i+1}(\xi) - A_i(\xi)]$ 方法2：加权最小二乘法（推荐）考虑到不同位置的采样质量不同，使用加权最小二乘： \[\Delta A_i = \arg\min_{\Delta} \sum_{\xi \in overlap} w(\xi) [A_{i+1}(\xi) - A_i(\xi) - \Delta]^2\] 其中权重 $w(\xi)$ 通常取为该点的采样次数：$w(\xi) = \min(N_i(\xi), N_{i+1}(\xi))$，确保采样好的区域有更高的权重。方法3：基于平均力的直接拼接更精确的方法是直接在重叠区域比较平均力，而非PMF： \[\Delta A_i = -\int_{\xi_{i+1}^{min}}^{\xi_i^{max}} [\langle F_{i+1}(\xi') \rangle - \langle F_i(\xi') \rangle] \mathrm{d}\xi'\] 这种方法对噪声更鲁棒，因为它利用了原始的平均力数据。步骤3：全局拼接从第一个窗口开始，逐步累积偏移量，构建全局PMF： \[A(\xi) = \begin{cases} A_1(\xi) & \text{if } \xi \in [\xi_1^{min}, \xi_1^{max}] \\ A_2(\xi) + \Delta A_1 & \text{if } \xi \in [\xi_2^{min}, \xi_2^{max}] \\ A_3(\xi) + \Delta A_1 + \Delta A_2 & \text{if } \xi \in [\xi_3^{min}, \xi_3^{max}] \\ \vdots \\ A_i(\xi) + \sum_{j=1}^{i-1} \Delta A_j & \text{if } \xi \in [\xi_i^{min}, \xi_i^{max}] \end{cases}\] 在重叠区域的处理：对于重叠区域 $[\xi_{i+1}^{min}, \xi_i^{max}]$，可以：选择其一：只使用窗口 $i$ 或窗口 $i+1$ 的数据加权平均（推荐）： $A(\xi) = \frac{w_i(\xi) \cdot [A_i(\xi) + \sum_{j=1}^{i-1}\Delta A_j] + w_{i+1}(\xi) \cdot [A_{i+1}(\xi) + \sum_{j=1}^{i}\Delta A_j]}{w_i(\xi) + w_{i+1}(\xi)}$ 其中 $w_i(\xi) = N_i(\xi)$ 是窗口 $i$ 在 $\xi$ 处的采样次数步骤4：质量检查拼接完成后，应检查：连续性：相邻窗口的PMF在重叠区域是否平滑连接一致性：重叠区域内两个窗口的PMF差异是否小于统计误差（通常 < 0.5 kcal/mol）平均力一致性：重叠区域内 $\langle F_i(\xi) \rangle$ 和 $\langle F_{i+1}(\xi) \rangle$ 是否接近与WHAM的对比： ABF拼接：简单、直接，只需在重叠区域对齐PMF，不需要迭代求解 WHAM：用于umbrella sampling，需要迭代求解自洽方程，计算复杂度更高，但在窗口重叠较少时更稳定四、ABF的优势与局限优势无需先验知识：不需要预先知道自由能曲面的形状高效采样：在能垒高的区域，ABF比umbrella sampling更高效无后处理：不需要WHAM等复杂的后处理方法局限初期采样问题：在模拟初期，平均力估计不准确，需要设置一个最小采样阈值（如每个网格点至少100次访问）才开始施加偏置隐藏能垒：如果正交于CV的自由度存在高能垒，ABF可能采样不充分几何约束的影响：当CV与几何约束或其他CV耦合时，需要使用扩展ABF（extended ABF, eABF）来正确处理五、主流MD软件中的ABF实现 5.1 NAMD中的ABF 实现方式：ABF在NAMD中通过Colvars模块（Collective Variables Module）实现，是NAMD内置的官方支持方法。基本使用流程：定义集合变量：在配置文件中定义CV（如距离、角度、二面角、RMSD等） colvar { name myDistance distance { group1 { atomNumbers 1 2 3 } group2 { atomNumbers 10 11 12 } } } 启用ABF：配置ABF参数 abf { colvars myDistance fullSamples 200 # 开始施加偏置前的最小采样数 historyfreq 50000 # 输出频率 writeTISamples yes # 输出统计数据 } 运行模拟：NAMD自动计算瞬时力、累积平均力并施加偏置支持的集合变量类型： distance：原子间距离 angle、dihedral：键角和二面角 rmsd：相对参考结构的RMSD gyration：回旋半径 eigenvector：沿主成分的投影输出文件： .pmf：PMF曲线数据 .count：每个网格点的采样次数 .grad：平均力数据参考资源： NAMD官方ABF教程：https://www.ks.uiuc.edu/Training/Tutorials/namd/ABF/ Colvars参考手册：https://colvars.github.io/colvars-refman-namd/ 5.2 GROMACS中的ABF 实现方式：GROMACS本身不直接支持ABF，但有以下几种替代方案：方案1：GROMACS + PLUMED（不推荐用于ABF） PLUMED是一个通用的增强采样插件，支持多种MD引擎局限：PLUMED不计算二阶导数，只能实现基于一阶导数的简化ABF版本 ABF并非PLUMED的原生方法，需要自行用C/C++实现方案2：GROMACS + SSAGES（推荐用于ABF） SSAGES（Software Suite for Advanced General Ensemble Simulations）提供了完整的ABF实现使用流程：使用GROMACS工具准备输入文件（拓扑、坐标）编写SSAGES的JSON配置文件定义CV和ABF参数使用gmx_ssages或gmx_mpi运行模拟文档：https://ssagesproject.github.io/ 方案3：GROMACS原生AWH方法（推荐替代） AWH（Accelerated Weight Histogram）是GROMACS 2018及以后版本的原生自适应偏置方法原理类似ABF：通过自适应调整偏置势来加速采样并计算PMF 优势： GROMACS原生支持，无需外部插件性能优化好，与GROMACS集成度高文档完善基本使用： pull = yes pull-ncoords = 1 pull-coord1-type = umbrella pull-coord1-geometry = distance pull-coord1-groups = 1 2 awh = yes awh-nstout = 1000 awh-nbias = 1 awh1-ndim = 1 awh1-dim1-coord-index = 1 参考文档：https://manual.gromacs.org/current/reference-manual/special/awh.html 推荐方案对比：方案优势劣势适用场景 SSAGES 完整ABF实现需要额外编译安装需要严格使用ABF算法 AWH 原生支持、性能好与标准ABF略有差异大多数自适应偏置应用 PLUMED 通用性强、功能多 ABF支持有限使用其他增强采样方法 5.3 其他MD软件 LAMMPS：通过Colvars模块支持ABF（与NAMD共用） Amber：通过PLUMED插件支持有限的ABF功能 OpenMM：通过Colvars或PLUMED插件支持总体建议：如需使用标准ABF方法，NAMD是首选（原生支持，文档完善） GROMACS用户建议使用AWH方法（原生、高效）或SSAGES（标准ABF）对于多维复杂CV或需要与其他增强采样方法结合，考虑使用PLUMED

Molecular Dynamics · 2025-10-09

跨越毫秒到秒级鸿沟：加权系综模拟如何捕捉”看不见”的生物动力学本文信息标题：加权系综模拟：方法、软件与应用的进展作者：Lillian T. Chong, Daniel M. Zuckerman 发表时间：2025年5月6日（ChemRxiv预印本）单位：匹兹堡大学（美国），俄勒冈健康与科学大学（美国）引用格式：Chong, L. T., & Zuckerman, D. M. (2025). WEIGHTED ENSEMBLE SIMULATION: ADVANCES IN METHODS, SOFTWARE, AND APPLICATIONS. ChemRxiv. https://doi.org/10.26434/chemrxiv-2025-jtppp 相关软件：本文主要讨论了基于 WESTPA 软件包的进展，并提及了其他实现如 wepy。摘要二十多年来，加权系综（Weighted Ensemble, WE）路径采样策略以远低于传统模拟的计算成本，实现了对罕见事件（或称跨能垒过程）路径的模拟，同时保持了严谨的动力学信息。本综述重点介绍了WE在方法和软件方面的最新进展，包括用于路径系综机理分析和高效速率估算的工具。我们展示了加权系综在一系列广泛的凝聚相过程中的成功应用，例如，微秒时间尺度的化学反应的混合量子力学/分子力学（QM/MM）模拟，以及毫秒到秒时间尺度的更慢过程的原子级模拟。这些应用涵盖了药物跨膜渗透、配体解离以及SARS-CoV-2刺突蛋白的大尺度开放等前沿领域。我们还讨论了WE策略当前面临的局限性和关键挑战，该方法尚未完全发挥其潜力。核心结论 WE是高效的罕见事件采样方法：它通过复制（分裂）和删减（合并）轨迹，能够以更低的计算成本模拟药物解离、蛋白质构象变化等低概率事件，同时严格保留动力学信息。方法学日趋成熟：近年来，WE在反应坐标优化（如机器学习辅助）、速率常数估算和不确定性量化等方面取得了显著进展，使其更加强大和可靠。软件生态系统完善：以 WESTPA 为代表的开源软件包具有高度可扩展性和互操作性，无需修改动力学引擎即可与AMBER、GROMACS、OpenMM等主流软件无缝对接，极大地促进了其应用。应用成果斐然：WE已成功应用于多个前沿领域，包括模拟秒级的SARS-CoV-2刺突蛋白开放、药物分子从深埋口袋中的解离、以及微秒级的QM/MM化学反应，揭示了实验难以企及的机理细节。背景 mindmap root((**背景与动机**)) **罕见事件采样挑战** 蛋白质折叠 **药物结合解离** 酶催化反应跨越能垒的过程 **传统MD模拟的限制** **毫秒到秒时间尺度** **计算成本高昂** 难以捕捉罕见事件 **WE策略的特色** **优胜劣汰重点培养** **动态资源分配** **无偏轨迹采样** 保留动力学信息在分子模拟的世界里，许多最引人入胜的生物学过程——如蛋白质折叠、药物分子与靶点的结合与解离、酶催化反应——都属于”罕见事件“。这意味着这些过程虽然至关重要，但在整个模拟时间尺度中，系统大部分时间都处于稳定的能量”盆地”中，而跨越能垒发生关键转变的瞬间则极为短暂和稀少。使用传统的分子动力学（MD）模拟，想要捕捉到这些事件的完整路径和动力学信息，往往需要运行长达毫秒、秒甚至更长时间的模拟，这对于目前的计算资源来说是极其昂贵甚至是不可能的。为了攻克这一难题，科学家们开发了多种增强采样和路径采样方法。其中，加权系综（Weighted Ensemble, WE）是一种尤为强大且独特的路径采样策略。与那些通过修改能量势面来加速转变的方法不同，WE的核心思想是”优胜劣汰，重点培养“。它并行地运行大量短时间的、完全标准的MD轨迹，并为每条轨迹分配一个”权重”。在固定的时间间隔后，它会评估所有轨迹的位置，智能地”克隆”那些正在向我们感兴趣的罕见区域探索的轨迹（分裂），并”删减”那些在已充分采样的区域中冗余的轨迹（合并）。通过这种方式，WE将计算资源动态地重新分配到那些”有前途”的路径上，极大地提高了采样到罕见事件的效率，同时由于每条轨迹本身是无偏的，整个过程保留了严谨的动力学信息，可以直接用来计算反应速率常数等关键物理量。经过二十多年的发展，WE方法本身、支持它的软件以及其应用范围都取得了长足的进步。关键科学问题作为一篇综述，本文旨在系统性地回答以下问题，为相关领域的研究者提供一份全面的指南和前沿展望： WE方法的核心原理是什么？它与其他路径采样方法相比有何独特的优势和固有的局限性？近年来WE方法学本身有哪些关键突破？研究者们是如何解决诸如如何定义”进展”、如何更准确地计算速率、以及如何评估结果不确定性等核心挑战的？支持WE模拟的软件生态系统发展如何？以WESTPA为代表的软件包在可扩展性、易用性和与其他主流模拟软件的兼容性方面取得了哪些进展？ WE在解决实际科学问题上取得了哪些里程碑式的应用成果？它如何帮助我们理解从病毒入侵到药物设计等一系列复杂生物过程的动力学机理？ WE方法的未来在哪里？它仍然面临哪些挑战，以及未来的发展方向将如何进一步拓展其应用边界？研究内容核心理论：加权系综（WE）模拟的”道”与”术” mindmap root((**WE核心原理**)) **基本算法** 初始化反应坐标定义箱子bins划分 **权重归一化** **动力学演化** 并行短时MD 无偏轨迹生成 **重采样操作** **分裂Splitting** 探索稀有区域克隆轨迹 **合并Merging** 删减冗余保持权重和为1 迭代循环 **动力学计算** **源-汇边界** **非平衡稳态NESS** **速率常数** **显著特点** 互操作性强算法灵活轨迹无偏连续统计严格精确 **固有局限** 物理时间尺度限制轨迹相关性问题方差挑战 WE方法的核心思想在于通过操控一个带有权重的轨迹系综，在不偏离真实动力学的前提下，高效地对罕见事件进行采样。基本算法流程：分裂与合并的艺术 WE方法通过并行管理一组带有权重的轨迹来高效探索构象空间。整个过程可以被看作一个不断迭代的循环。图1：加权系综策略示意图该图展示了一个基础的WE实现，其中构象空间被划分为固定的”箱子（bins）”，每个箱子的目标轨迹数为3。子轨迹会均分其父轨迹的权重，确保每一轮迭代中总权重为1。初始化与空间划分：首先，需要定义一个或多个”反应坐标（Progress Coordinates）“，它们是能够描述系统从初始态向目标态演化进程的变量。基于这些坐标，整个构象空间被划分成一系列离散的”箱子（bins）“。然后，从一个或多个初始构象开始，启动若干条轨迹，并为它们分配初始权重。所有轨迹的权重总和必须恒为1，即： \[\sum_{i} w_i(t) = 1\] 动力学演化（Evolve）：在一个迭代步中，所有轨迹都独立、无偏地进行一小段固定时间（$\tau$）的MD模拟。这个步骤是完全并行的，因此WE具有极好的可扩展性。重采样（Resampling）：这是WE的灵魂所在。在 $\tau$ 时间后，暂停所有轨迹，并根据它们所处的”箱子”进行分裂（Splitting）和合并（Merging）操作：分裂（Splitting）：当一条轨迹进入了一个很少被访问或完全空的箱子时，表明它正在探索新的、重要的区域。此时，系统会将其”克隆”成两条或多条（例如2条）子轨迹。这些子轨迹完全继承父轨迹的坐标和速度，并均分其权重（例如，权重为 $w_p$ 的父轨迹分裂成两条权重各为 $w_p/2$ 的子轨迹）。这相当于将计算资源动态地聚焦到有前途的探索路径上。合并（Merging）：当一个箱子里的轨迹数量超过了预设的目标值时，说明该区域已被过度采样，存在冗余计算。此时，系统会从中选择轨迹进行合并。例如，从箱子中随机选取两条轨迹 $i$ 和 $j$，它们的权重分别为 $w_i$ 和 $w_j$。系统会根据权重以概率 $p_i = w_i / (w_i + w_j)$ 保留轨迹 $i$，或以概率 $p_j = w_j / (w_i + w_j)$ 保留轨迹 $j$。幸存的轨迹将获得两者合并后的总权重 $w_{\text{new}} = w_i + w_j$，而被淘汰的轨迹则终止。这相当于剪除冗余的计算分支，节约资源。迭代：完成重采样后，所有”幸存”和”新生”的轨迹进入下一轮的”演化-重采样”循环，周而复始，直到达到预定的总模拟时间或目标事件被充分采样。 graph TD subgraph "方向：从左到右" direction LR A("1.初始化 一组带权重的轨迹") --> B["2.动力学演化 所有轨迹独立运行一小段时间 τ"]; B --> C{"3.重采样 (根据轨迹位置)"}; C -- "进入稀有区域" --> D["分裂 (复制有前途的轨迹)"]; C -- "进入拥挤区域" --> E["合并 (删除冗余的轨迹)"]; D --> F["进入下一轮迭代"]; E --> F; end 动力学性质计算：速率常数 WE的一个核心优势是能够直接计算动力学速率常数。这通常通过设置”源-汇（source-sink）“边界条件来实现：当一条轨迹到达我们定义的目标态（汇），它不会终止，而是被”传送”回初始态（源）并继续模拟。经过一段时间的模拟，系统会达到一个非平衡稳态（Non-Equilibrium Steady State, NESS），此时单位时间内从源到达汇的概率通量（Flux）将趋于一个稳定值，这个值就是我们要求的速率常数 $k_{AB}$。 \[k_{AB} = \text{Flux}(A \rightarrow B | \text{NESS})\] 公式的通俗解释这个公式是WE计算速率的核心。 $k_{AB}$：是从状态A到状态B的速率常数，单位是时间的倒数（如 $\mathrm{s}^{-1}$）。 $\text{Flux}(A \rightarrow B)$：指的是单位时间内，从初始态A区域”流向”目标态B区域的总概率。在WE中，这就是所有首次到达目标态B的轨迹的权重之和除以时间间隔 $\tau$。 NESS：表示这个计算必须在系统达到非平衡稳态后进行。如图2所示，模拟刚开始时，通量会逐渐增加（瞬态），只有当进入和离开各个区域的概率流达到一种动态平衡时，测得的通量才是稳定且准确的。图2：从WE模拟流入目标态的通量估计速率常数模拟开始后，流入目标态的通量会经历一个瞬态增长期，最终达到一个平台期，即非平衡稳态，此时的通量值即为速率常数 $k$。 WE方法的显著特点与固有局限优点互操作性强：WE算法只要求能启停轨迹，因此无需修改任何MD引擎的底层代码，可以与AMBER、GROMACS、OpenMM等任何模拟软件无缝协作。这种设计使得研究者可以继续使用最熟悉、最适合其体系的MD引擎，而不必为了使用WE而去学习一个全新的模拟软件。算法灵活：WE的分箱策略、资源分配等都可以在模拟过程中动态调整，甚至可以完全抛弃”箱子”概念，而是基于轨迹间的相似度进行重采样（如REVO方案）。这种高度的灵活性使得WE能够适应各种复杂的生物分子体系。轨迹无偏且连续：WE不施加任何偏置力，每条轨迹片段都是真实的动力学路径，最终可以拼接成完整的、可用于各种机理分析的连续轨迹。这种无偏性是WE与其他增强采样方法的根本区别。统计上严格精确：理论上，WE的系综平均结果与大量传统MD模拟的结果是完全一致的。这种统计上的严谨性使得WE计算得到的速率常数等动力学可观测量具有理论上的精确性。高效并行性：WE具有极好的可扩展性，能够在数千个CPU/GPU核心上高效并行，其优秀的任务管理器和通信层设计使其能够驾驭超算级别的计算资源。计算成本显著降低：相比传统MD模拟，WE能够以远低于传统模拟的计算成本实现对罕见事件（或称跨能垒过程）路径的模拟，同时保持严谨的动力学信息。局限性 2.3 WE的内在局限性 WE方法的主要局限性源于系统固有的物理时间尺度，因此这是任何模拟真实连续轨迹系综的方法都面临的挑战。具体而言，任何感兴趣的转变过程都可以用平均过渡路径时间（average transition path time） $\langle t_{\text{TP}} \rangle$ 来表征。因此，包含 $n \gg 1$ 条轨迹的系综所需的总时间为 $n \cdot \langle t_{\text{TP}} \rangle$，这代表了在能够完全独立生成正确分布的过渡轨迹（这实际上是不可能的）的理想情况下的绝对最小计算成本。实际上，还存在一个额外的低效因子 $m > 1$（很可能 $m \gg 1$），它代表了生成独立轨迹的开销成本。因此，系综的总成本为 $m \cdot n \cdot \langle t_{\text{TP}} \rangle$，这还没有考虑 $t_{\text{TP}}$ 值可能存在的非高斯大幅度涨落。即使对于 $\langle t_{\text{TP}} \rangle \sim 10 \mathrm{ns}$ 的转变过程，在WE或其他生成连续路径系综的方法中也可能需要数微秒的轨迹数据。对 $\langle t_{\text{TP}} \rangle$ 的估计各不相同：小蛋白折叠（微秒到毫秒时间尺度）约为1-100 ns，扩散控制的蛋白-蛋白结合（微秒时间尺度）约为5 ns，蛋白-配体解离（秒时间尺度）约为100 ns。为什么高度相关轨迹会导致WE估计的可观测量（如速率常数）在不同运行之间存在高方差？统计独立性缺失：在WE中，分裂操作产生的子轨迹共享相同的历史，导致它们高度相关。这些相关轨迹不提供独立的统计信息，相当于减少了有效样本量。当多个相关轨迹贡献到同一统计量时，它们不能像独立轨迹那样有效降低方差，导致估计的不确定性增加。路径空间采样不均衡：相关轨迹倾向于探索相似的路径空间区域，使得某些重要但罕见的路径可能被低估，而常见路径则被过度采样。这种采样不均衡性会导致不同WE运行之间对同一物理量的估计出现较大波动。权重分布偏差：由于合并操作基于权重进行随机选择，高度相关的轨迹可能导致权重分布出现偏差。这种权重偏差会进一步放大估计量的方差，尤其是在长时间模拟中。收敛速度降低：相关轨迹减慢统计收敛速度，因为系统需要更长时间探索不同的路径空间。在有限的计算资源下，这可能导致不同运行之间结果差异显著。低效因子 $m$ 正是反映了这种基于相关性的低效率。在WE中，这种相关性源自基本的分裂操作。同一父轨迹的子轨迹在分裂点之前共享相同的历史，使它们高度相关。高度相关轨迹的实际后果是WE估计的可观测量（如速率常数）在不同运行之间可能存在高方差（图2）。这种基于相关性的方差和低效率可以在一定程度上得到改善，下文将详细讨论。我们还注意到，相关性使得不确定性量化更具挑战性，这也将在下文讨论。总体而言，虽然WE是一种强大而严格的方法，但并不保证在每个系统上都能很好地工作。例如，高电荷配体从蛋白受体解离是一个特别具有挑战性的压力测试；相比之下，更容易处理的应用涉及不带电配体的解离（见第5.2节）。基于系统物理性质的固有成本是显著的，这不仅对WE如此，对任何提供真实过渡路径系综的方法都是如此，即使是粗粒化模型也是如此。基于相关性的低效率也是路径采样方法的固有特性。 WE方法学的最新进展图3：WE方法学中的挑战与解决方案（a）WE模拟面临着寻找好的反应坐标、速率估计方差大和不确定性量化等挑战。（b）针对这些挑战，研究者开发了机器学习、方差最小化分箱和贝叶斯分析等解决方案。（c）这些优化方法通常需要初步模拟数据，通过分析或机器学习来指导后续的优化模拟或直接计算可观测量。近年来，研究者们从多个角度对WE方法进行了优化，主要分为两大类：优化模拟过程：反应坐标与分箱策略：这是WE实践中最关键的一环。除了依赖化学直觉，多种自动化策略被开发出来。例如，最小自适应分箱（MAB）方案能自动识别路径上的瓶颈区域并增加采样；REVO 方案则完全抛弃箱子，基于轨迹相似度进行重采样。机器学习也被用于从业已产生的高维轨迹数据中自动学习出最优的低维反应坐标，例如使用卷积变分自编码器来压缩轨迹信息。更有甚者，可以直接以最小化速率常数估计的方差为目标来优化分箱策略。优化数据分析：速率常数估算：为了解决模拟时间不足以达到稳态的问题，研究者开发了历史增强马尔可夫状态模型（haMSM），它可以从非稳态的瞬态数据中外推出稳态的速率常数。机理量化：如何从大量的路径中提取并量化”反应机理”是一个开放性问题。目前已有如LPATH等工具被开发用于对路径进行聚类和分析，以识别不同的反应通道。不确定性量化（UQ）：由于轨迹相关性，简单的统计方法不适用。目前的主流做法是进行多次独立的WE模拟，然后分析多次模拟结果之间的差异，有时还会借助贝叶斯分析来处理方差较大的情况。 WE软件的进展：以WESTPA为例 mindmap root((**WESTPA软件生态**)) **可扩展性** 数千CPU和GPU核心超算级别支持优秀任务管理器通信层设计 **互操作性** **与引擎解耦** 命令行调用 **支持主流软件** AMBER GROMACS OpenMM 无需代码修改 **数据管理** **WESTPA2.0改进** **HDF5格式** 高效存储便利重启分析数据共享优化 **未来发展** Dask任务分发减少延迟容错能力云计算支持 WESTPA（The Weighted Ensemble Simulation Toolkit with Parallelization and Analysis）是目前最活跃、功能最强大的开源WE软件包之一。高度可扩展：WESTPA能够在数千个CPU/GPU核心上高效并行，其优秀的任务管理器和通信层设计使其能够驾驭超算级别的计算资源。强大的互操作性：WESTPA设计上与动力学引擎解耦，可以像”指挥官”一样通过命令行调用任何模拟软件（如AMBER、GROMACS、OpenMM）或分析工具（如MDAnalysis、MDTraj），无需任何代码修改。数据管理优化：最新的WESTPA 2.0版本改进了数据存储框架，使用高效的HDF5格式来管理数千万个轨迹文件，极大地便利了模拟重启、数据共享和后分析。未来发展：未来的WESTPA将集成更先进的任务分发框架（如Dask），以减少延迟、增强容错能力，并更好地支持云计算平台。 WE应用的亮点成果 mindmap root((**WE应用领域**)) **病毒学** **SARSCOV2刺突蛋白** 秒级时间尺度百万原子体系 **戈登贝尔奖** 聚糖门控机制实验验证 **药物发现** **配体解离** 秒级过程不带电配体 **HIF2α靶点** 两条解离路径 **隐蔽口袋探索** 不可成药靶点药物设计新路线 **跨膜渗透** 虚拟生物利用度 **ADMETOX评估** POPC脂双层与实验一致机理洞察 **化学反应** **QMMM模拟** 微秒级反应 click反应颠覆扩散控制假设限速步骤分析 **蛋白质相互作用** 结合动力学 kon速率计算关键残基识别蛋白质折叠图4：近期WE在微秒至秒时间尺度上的应用（a）微秒级：化学反应的QM/MM模拟。（b）毫秒级：药物分子的跨膜渗透。（c）秒级：配体从深埋的受体口袋中解离。（d）秒级（百万原子体系）：SARS-CoV-2刺突蛋白的开放过程。病毒学：SARS-CoV-2刺突蛋白开放迄今为止最雄心勃勃的WE应用是对包含近百万个原子的SARS-CoV-2刺突蛋白（S蛋白）开放过程的模拟，这是一个秒级时间尺度的事件。这项工作荣获了2020年戈登·贝尔COVID-19研究特别奖。模拟不仅捕捉到了S蛋白从”关闭”到”开放”状态的完整路径，还揭示了一个前所未知的机理：位于N288位点的一个聚糖扮演了”构象门”的角色，控制着蛋白的开放。这一发现随后得到了实验的验证，包括生物层干涉测量实验和冷冻电镜（采用ManifoldEM方法生成S蛋白的大尺度运动，发现与模拟一致）。药物发现：配体解离与”隐蔽口袋”探索药物的疗效与其在靶点上的停留时间（与解离速率成反比）密切相关。WE已被成功用于模拟药物分子从靶点蛋白解离的秒级过程，迄今限于不带电配体。在一项针对癌症靶点HIF-2α PAS-B结构域的研究中，WE模拟在不知道任何先验信息的情况下，成功捕捉到了一个药物样不带电配体从其深埋的内部口袋中逃逸的两条不同路径。这些路径是以盲目方式生成的，无需任何关于解离过程的先验知识。模拟发现的构象门控残基也得到了NMR动力学实验的证实。此外，WE还能采样到在实验结构中不可见的”隐蔽口袋”，为”不可成药”靶点提供了潜在的可行药物设计路线。药物跨膜渗透：虚拟生物利用度分析 WE被用于开发预测药物被动跨膜渗透性的”虚拟实验”，这是评估药物吸收、分布、代谢、排泄和毒性（ADME/Tox）的关键性质。作为概念验证，WE模拟评估了一系列不同大小、形状和柔性的药物样胺类化合物通过模型POPC脂双层的渗透性。结果产生的渗透系数与MDCK-LE细胞系和平行人工膜渗透实验（PAMPA）的实验值一致，同时提供了转运过程的机理洞察。值得注意的是，尽管使用了被其他方法认为次优的反应坐标（膜中的z位置），WE仍成功生成了路径和速率估计，计算成本比传统MD低几个数量级。因此，WE策略对反应坐标选择的敏感性远低于基于自由能的方法。化学反应：QM/MM模拟揭示反应机理通过与混合量子力学/分子力学（QM/MM）方法结合，WE首次被用于模拟溶液中的化学反应并计算速率。在一项对叠氮化物”click反应”的研究中（叠氮阴离子与三苯甲基阳离子在乙腈-水溶液中反应），WE-QM/MM模拟不仅重现了实验速率，还颠覆了之前的”扩散控制”假设，指出反应的限速步骤是离子对中间体重排为产物的活化过程。研究还揭示了叠氮离子在阳离子苯环间增加的”爬行”与更慢的反应速率相关，这项工作突显了WE在使用混合QM/MM模型进行路径采样和动力学分析以获得更深入机理洞察方面的威力。蛋白质-蛋白质相互作用：结合动力学研究 WE已被用于研究蛋白质-蛋白质结合路径和速率常数计算，采用完全连续的显式溶剂模拟。通过模拟结合路径和解离过程，WE能够揭示相互作用界面的关键残基和构象变化。例如，WE已被用于计算基础 $k_{\text{on}}$（直接模拟柔性分子模型的蛋白-蛋白结合），以及比较无序肽及其精确预组织类似物的结合动力学。蛋白质折叠：超快折叠蛋白研究 WE已被成功应用于研究蛋白质折叠动力学和机制。例如，在对超快折叠蛋白NTL9的研究中，WE模拟揭示了改变骨架组成对折叠动力学和机制的影响。这些应用展示了WE在解决从微秒到秒原子级折叠时间的计算估计方面的能力。多尺度过程与未来展望 WE方法的应用范围正在不断扩展。除了上述应用，WE还被应用于肽跨膜渗透、脂质相分离热力学、以及大规模生物分子复合物的动力学研究。随着计算能力的提升和方法的持续改进，WE有望在更复杂的细胞环境（如呼吸道气溶胶、细菌或人类细胞质）中模拟生物分子的行为。 Q&A Q1：加权系综（WE）和其他增强采样方法（如元动力学、伞形采样）的根本区别是什么？ A1：根本区别在于是否改变系统的哈密顿量（即能量势面）。元动力学、伞形采样等方法属于偏置势（Biasing Potential）方法。它们通过在构象空间中添加一个外部的、人为的偏置势能来”填平”能量势垒，从而迫使系统更快地在不同状态间转换。这些方法能高效地计算自由能曲线，但其产生的轨迹不是真实的动力学路径，因此不能直接用来计算速率常数或分析动力学机理。加权系综（WE）则是一种路径采样（Path Sampling）方法。它不施加任何偏置力，系统在每一步都遵循自然的动力学演化。它的加速效果来自于在路径空间中对轨迹进行智能的复制和删减，即把计算资源集中到更有可能发生转变的路径上。因此，WE产生的轨迹是物理上真实的、无偏的连续路径，既可以用来计算自由能，也可以直接用来分析动力学机理和计算速率常数。 Q2：什么是好的”反应坐标（progress coordinate）”，为什么它对WE模拟如此重要？ A2：一个好的”反应坐标”是一个或一组能够有效区分反应物、产物以及过渡态的低维变量。它应该能够捕捉到系统从初始态向目标态的”进展程度”。在WE模拟中，反应坐标直接决定了”箱子（bins）”的划分，从而控制着轨迹的分裂与合并策略。一个好的反应坐标能让WE算法准确地识别出哪些轨迹正在接近反应的”瓶颈”区域（即能垒顶部），并及时在这些关键区域增加采样（分裂轨迹），从而大大提高模拟效率。相反，如果选择了一个与反应真实路径无关的坐标，WE可能会在不相关的区域浪费大量计算资源，导致收敛缓慢甚至失败。 Q3：WESTPA软件的一大亮点是”互操作性（interoperability）”，这具体指什么，为什么它很重要？ A3：互操作性指的是WESTPA能够与几乎任何现有的动力学模拟软件（如AMBER、GROMACS、OpenMM）或分析工具无缝协作，而无需对这些软件进行任何代码修改。WESTPA就像一个”总指挥”，它通过标准的命令行接口来启动、监控和停止由其他软件执行的短时间模拟任务，然后在每个迭代周期结束后收集结果并执行重采样。这一点至关重要，因为它极大地降低了使用WE方法的门槛。研究者可以继续使用他们最熟悉、最适合其体系的MD引擎，而不必为了使用WE而去学习一个全新的、功能可能不全的模拟软件。这种模块化的设计也使得更换动力学引擎或升级版本变得非常简单。关键结论与批判性总结潜在影响解锁长时程动力学：WE及其相关软件的发展，使得在原子级别上直接模拟并分析毫秒至秒级甚至更长时间尺度的生物过程成为可能，为理解药物停留时间、病毒入侵机理等关键问题提供了前所未有的工具。连接理论与实验：WE能够直接计算速率常数等动力学可观测量，这为力场的动力学性质验证提供了黄金标准，有助于推动下一代更精确的分子力场的开发。推动多尺度模拟：WE的灵活性使其不仅限于分子模拟，还可以应用于系统生物学、天气预报等更宏观的尺度，展现了其作为一种通用罕见事件采样方法的巨大潜力。研究局限性方法仍在发展中：尽管取得了巨大成功，但WE方法仍处于活跃的发展阶段。如何系统性地选择最优反应坐标、如何更精确地进行不确定性量化等问题仍是当前研究的热点和挑战。对特定体系的挑战：对于某些体系，如高电荷配体的解离，WE模拟仍然面临巨大挑战，结果的方差可能非常大，难以收敛。资源需求依然可观：虽然WE相比传统MD效率极高，但模拟秒级过程仍然需要巨大的计算资源（如SARS-CoV-2的研究），这限制了其在普通实验室的广泛应用。未来方向 QM/MM与WE的深度融合：进一步推动WE在QM/MM模拟中的应用，有望在更长的时间尺度（多微秒级）上研究酶催化和溶液中的化学反应。超长时程模拟：随着计算能力的提升和算法的持续优化，WE有望挑战秒级以上的生物过程，为研究治疗性相关的动力学事件提供更精确的速率估计。与实验数据的整合：将WE产生的路径系综与单分子实验（如FRET）或时间分辨结构生物学数据相结合，以更全面的视角揭示生物大分子的功能机理。向更复杂环境迈进：随着细胞环境的结构数据日益丰富，未来的WE模拟将不再局限于孤立的生物分子，而是能够模拟其在呼吸道气溶胶、细菌乃至人类细胞质等更真实、更拥挤环境中的行为。

Molecular Dynamics

Mendelevium

Contact

Fundamental